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Abstract 

We have formulated the problem of generating dense packings of nonoverlapping, non-tiling 
nonspherical particles within an adaptive fundamental cell subject to periodic boundary conditions 
as an optimization problem called the Adaptive Shrinking Cell (ASC) formulation [S. Torquato 
and Y. Jiao, Phys. Rev. E 80, 041104 (2009)]. Because the objective function and impenetrability 
constraints can be exactly linearized for sphere packings with a size distribution in (i-dimensional 
Euclidean space R rf , it is most suitable and natural to solve the corresponding ASC optimization 
problem using sequential linear programming (SLP) techniques. We implement an SLP solution 
to produce robustly a wide spectrum of jammed sphere packings in W 1 for d = 2, 3, 4, 5 and 6 with 
a diversity of disorder and densities up to the respective maximal densities. A novel feature of 
this deterministic algorithm is that it can produce a broad range of inherent structures (locally 
maximally dense and mechanically stable packings), besides the usual disordered ones (such as the 
maximally random jammed state), with very small computational cost compared to that of the 
best known packing algorithms by tuning the radius of the influence sphere. For example, in three 
dimensions, we show that it can produce with high probability a variety of strictly jammed packings 
with a packing density anywhere in the wide range [0.6,0.7408...], where 7r/\/l8 = 0.7408... 
corresponds to the density of the densest packing. We also apply the algorithm to generate various 
disordered packings as well as the maximally dense packings for d = 2,4, 5 and 6. Our jammed 
sphere packings are characterized and compared to the corresponding packings generated by the 
well-known Lubachevsky-Stillinger molecular-dynamics packing algorithm. Compared to the LS 
procedure, our SLP protocol is able to ensure that the final packings are truly jammed, produces 
disordered jammed packings with anomalously low densities, and is appreciably more robust and 
computationally faster at generating maximally dense packings, especially as the space dimension 
increases. 

PACS numbers: 05.20.Jj, 45.70.-n, 61.50.Ah 
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I. INTRODUCTION 



Hard-particle packings have provided a rich source of outstanding theoretical problems 



and served as useful starting points to model t 



Q, QJ18 



crystals 7H9|, colloids [2|, |10J, granular media 



re structure of liquids 
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and polymeric systems 
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14J, living cells 



l[ 3, glasses |a-fj], 
1511. random media 



20|. We focus our attention in this paper on the 
venerable idealized hard-sphere model in (i-dimensional Euclidean space M. d in which the 
only interparticle interaction is an infinite repulsion for overlapping particles, which can be 
thought of as the "Ising model" for hard spheres 

There has been resurgent interest in hard-sphere packings in dimensions greater than 
three in both the physical and mathematical sciences. For example, it is known that the 
optimal way of sending digital signals over noisy channels corresponds to the densest sphere 
packing in a high-dimensional space 
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221 ] . These error-correcting codes underlie a variety 



of systems in digital communications and storage, including compact disks, cell phones 
and the Internet. Physicists have studied hard-sphere packings in high dimensions to gain 
insight into ground and glassy states of matter as well as phase behavior in lower dimensions 
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. The determination of the densest packings in arbitrary dimension is a problem of 
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3Q|. 



longstanding interest in discrete geometry and number theory 

The hard-sphere Ising model enables us to be precise about the important concept of 
"jamming." The generation and characterization of jammed packings are the main concerns 
of this paper. A jammed packing is one in which the particle positions are fixed by the 
impenetrability constraints and boundary conditions (e.g., hard- wall or periodic boundary 
conditions). For sphere packings, one can define hierarchical jamming categories beginning 



from the least restrictive to the most stringent one [13|. In particular, a packing is locally 
jammed if no particle in the system can be translated while fixing the positions of all other 
particles, requiring that each sphere in the packing be in contact with at least d + 1 spheres 
not all in the same hemisphere. A collectively jammed packing is a locally jammed packing 
such that no subset of spheres can simultaneously be continuously displaced so that its 
members move out of contact with one another and with the remainder set. A packing is 
strictly jammed if it is collectively jammed and all globally uniform volume nonincreasing 
deformations of the system boundar y a re disallowed by the impenetrability constraints. The 
reader is referred to Refs. 13| and 14J for further details. Our main concern in this paper 
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will be packings that are at least collectively jammed. 

The packing density </>, the fraction of M. d covered by the spheres, is the simplest char- 
acteristic of a jammed packing. However, such a characterization is clearly not sufficient 
in order to distinguish between ordered and disordered packings. In fact, jammed packings 
may be produced with variable degrees of disorder /order. It has been suggested that a 
scalar order metric ip be employed to measure the degree of order in a packing, such that 
ip = 1 corresponds to fully ordered [e.g., the perfect face-centered cubic (fee) crystal in 
three dimensions] and ip = corresponds to perfectly a disordered Poisson distribution of 



sphere centers 



311 ] . Thus, very large jammed packings have been classified based on their 



locations in the density-order (p-ip plane, called the order map. The points on the boundary 
of the jammed domain in the order map constitute optimal packings. For example, in M 3 , 
the maximally dense fee lattice sphere packing ((f) = 7r/vT8 = 0.7048 . . .) 32[ also has the 



highest order metric (ip = I) when appropriately defined. Another extremal state of special 
interest is the maximally random jammed (MR J) state, which is the one that minimizes a 



scalar order metric ip s ubj ect to the condition of the degree of jamming 



Studies of 



different order metrics 33| for three-dimensional frictionless spheres have consistently led to 
a minimum at approximately the same density <ft ~ 0.64 for collective and strict jamming 
in the order map (4, 31] . Such MRJ packings have also been shown to be hyperuniform, 

34j. The lowest density strictly jammed 



i.e., infinite- wave-length density fluctuations vanish 

states in M 3 , thought to be tunneled crystals js], are a fascinating set of extremal loci in 
the jamming region of the order map [14]. The frequency of occurrence of a particular 
configuration is irrelevant insofar as the order map is concerned, i.e., the order map em- 
phasizes a "geometric-structure" approach to analyze packings by characterizing individual 
configurations, regardless of their occurrence probability l^ . 

During the last two decades, the Lubachevsky-Stillinger (LS) algorithm 35( has been 
the premier workhorse to generate a wide spectrum of dense jammed sphere packings with 



variable disorder in both two and three dimensions |4J, [311, |33|. This is an event-driven (or 
collision-driven) molecular dynamics algorithm in which an initial configuration of spheres 
of a given size within a periodic box are given initial random velocities and the motion 
of the spheres are followed as they collide elastically and also expand uniformly until the 
spheres can no longer expand. This algorithm has been generalized by Donev, Torquato 
and Stillinger 36[ to generate jammed packings of smoothly-shaped nonspherical particles, 
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including ellipsoids [37|, superdisks [38| and superballs 

Not surprisingly, this packing protocol is not without some inadequacies. Event-driven 
packing protocols with growing particles do not guarantee jamming of the final packing 
configuration, since jamming is not explicitly incorporated as a termination criterion. To 
produce jammed random packings, for example, a large expansion rate is necessary in the 
early stages of the simulation to suppress crystallization, but using a high expansion rate 
is highly undesirable toward the end of the simulation, which often leads to unjammed 
configurations. Thus either a variable expansion rate must be used (which decreases as a 
function time in some arbitrary fashion) or, if a uniform expansion is employed, the spheres 
of the terminal packing must be shrunk by some arbitrary small amount and then this initial 
packing configuration must be re-densified using a very small expansion rate. Maximally 
dense jammed packings are highly computationally expensive to generate, especially in high 
dimensions, because a very small expansion rate is necessary (on the order of 10" 6 to 10" 10 



401 ] . dependent on dimension). We will see that in five and six dimensions, even a expansion 
rate of 1CT 10 fails to produce the maximal-density packings. Moreover, a large number of 
total collisions per particle is required; for MRJ packings and the densest known maximally- 
dense packings, on the order of 10 5 and 10 7 collisions per particle are required, respectively, 
the latter of which is computationally very costly. Finally, such event-driven packing proto- 
cols evolve stochastic velocity initializations, which makes one have less control of the final 
packings via the initial configuration. 

The next generation packing protocol to generate jammed sphere packings should re- 
tain the versatility and advantages of the LS algorithm while correcting its imperfections, 
including improving computational speed. We show that our proposed sequential linear 
programming (SLP) solution of the Adaptive Shrinking Cell (ASC) optimization problem 
formulated elsewhere for general particle shapes (including polyhedra) [43( indeed has all 
of these desirable features in so far as jammed sphere packings with a size distribution are 
concerned. Because the design variables (including periodic simulation box shape and size) 
and impenetrability constraints can be exactly linearized for spheres packings, the deter- 
ministic SLP solution in principle always leads to jammed packings (up to a high numerical 
tolerance) with a wider range of densities and degree of disorder than packings produced 
by the LS algorithm. Each linear programming solution step starting from some initial 
particle configuration involves a deterministic collective motion of the entire particle con- 
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figuration to a higher density and, because the periodic simulation box can simultaneously 
deform and shrink, the final state is guaranteed to be at least collectively jammed in princi- 
ple. Whereas the LS algorithm requires between 10 5 and 10 7 collisions per particle without 
ensuring jammed final states, the deterministic SLP solution, which is easy to implement, 
requires only 10 to 100 linear-programming steps to achieve jammed packings, depending 
on the desired density, which can be controlled by tuning what we call the size of the in- 
fluence sphere. Thus, the SLP algorithm is computationally very efficient in generating the 
maximally dense packings, even in high space dimensions. By appropriately choosing initial 
conditions, one can achieve strictly jammed disordered sphere packings with anomalously 
low densities, e.g. <p ~ 0.6 in three dimensions 4l|, |42|. Indeed, a novel feature of the algo- 



rithm is that it can produce a broad range of inherent structures (locally maximally dense 
and mechanically stable packings), such as the maximally random jammed states as well as 
the globally maximally dense inherent structures, with very small computational cost. 

The rest of the paper is organized as follows. In Sec. II, we provide basic definitions 
for packing problems. In Sec. Ill, we present the mathematical formulation and algorith- 
mic details of our sequential-linear-programming (SLP) procedure to solve the adaptive- 
shrinking-cell (ASC) optimization problem to generate jammed sphere packings. In Sec. IV, 
we discuss the energy landscape (negative of the packing density) for jammed sphere pack- 
ings (inherent structures) and show that our SLP solution procedure is able to lead to both 
the mechanically stable local minima (i.e., local density maxima) and the global minima 
(maximal-density packings). In Sec. V, we employ the SLP algorithm to produce a diverse 
class of disordered jammed packings, including MRJ states, as well as maximal-density 
packings of hard spheres in M. d for d = 2, 3, 4, 5 and 6. The characteristics of these packings 
are compared to those generated using the LS algorithm. In Sec. VI, we employ the SLP 
algorithm to produce jammed packings in R 3 with varying degrees of disorder and packing 
densities anywhere in the wide range [0.6, 0.7408 . . .]. In Sec. VII, we discuss the ramification 
of our results and make concluding remarks. 

II. DEFINITIONS 

A lattice A in M d is a subgroup consisting of integer linear combinations of vectors Aj 
that constitute a basis for R^. In a lattice, the space M d can be geometrically divided into 
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identical regions F called fundamental cells, each of which contains one lattice site specified 
by the lattice vector 

p = mAi + n 2 X 2 H hn d _iA d _i +n d X d , (1) 

where Aj (i = 1, 2, . . . , d) are the basis vectors and n« spans all the integers for i — 1, 2, • • • d. 
The generator matrix Ma = {Ai, A2, . . . , A^} of a lattice A is a matrix with the basis vectors 
Aj as columns and, thus, contains d 2 elements. In the physical sciences and engineering, a 
lattice A is referred to as a Bravais lattice. Unless otherwise stated, the term "lattice" will 
refer here to a Bravais lattice only. A lattice packing of congruent spheres Pl is one in which 
the centers of the nonoverlapping identical spheres are located at the points of A and hence 
each fundamental cell contains within it exactly one sphere. Thus, the density of a lattice 
packing of spheres of radius R is given by 

*=^, (2) 

v F 

where vp is the volume of the fundamental cell and 

7T d ' 2 R" 

= F(ITW (3) 

is the volume of a d- dimensional spherical particle of radius R and T(x) is the Euler gamma 
function. 

A more general notion than a lattice packing is a periodic packing of spheres (not nec- 
essarily identical in size), which is obtained by placing a fixed configuration of N spheres 
(where N > 1) of radii Ri, R 2 , ■ ■ ■ , Rn in one fundamental cell of a lattice A, which is then 
periodically replicated without overlaps. Thus, the packing is still periodic under transla- 
tions by A, but the N spheres can occur anywhere in the chosen fundamental cell subject 
to the overall nonoverlap condition. The packing density of a periodic packing of spheres is 
given by 

(j)= El 1 v 1 (R l ) = {vi{R)) (4) 
v F 

where p = N/vp is the number density, i.e., the number of particles per unit volume, and 

1 N 

(MR)) = jjJ2 v iW ( 5 ) 

i=l 

expected volume of a sphere. Figure Q] illustrates in IR 2 a lattice packing of spheres (which 
by definition requires that they be identical is size) and a periodic packing of spheres in 
which the spheres have a size distribution. 



(a) (b) 



FIG. 1: (Color online). Illustrations of two basic sphere packings in M?: (a) lattice packing of 
circular hard disks in M 2 , which necessarily can only involve congruent circles, and (b) periodic 
packing of circles of different sizes. 

III. MATHEMATICAL AND ALGORITHMIC DETAILS 

We have formulated the problem of generating dense packings of nonoverlapping, non- 
spherical particles within a fundamental cell in M 3 subject to periodic boundary conditions 



as an optimization problem called the the Adaptive Shrinking Cell (ASC) scheme 43j. The 
objective function is taken to be the negative of the packing density <ft. Starting from an 
initial unsaturated packing configuration of particles of fixed size in the fundamental cell, the 
positions and orientations of the particles are design variables for the optimization. We also 
allow the boundary of the fundamental cell to deform macroscopically as well as compress 
or expand (while keeping the particles fixed in size) such that there is a net compression 
(increase of the density of the packing) in the final state. Therefore, the deformation and 
compression/expansion of the cell boundary, called the adaptive fundamental cell, are also 
design variables. 

The ASC optimization problem can be solved using various techniques, depending on 
the shapes of the particles. For example, for polyhedron particles the nonoverlap condition 
is highly nonlinear, which makes it inefficient to solve even using nonlinear programming 
methods, and hence Monte Carlo techniques provide an efficient means of solving the problem 



8 



[43]. However, for spheres, one can exactly linearize the objective function, design variables 
and constraints, enabling one to exploit efficient linear programming techniques. 

We first explicitly formulate the ASC scheme for noncongruent sphere packings, i.e., 
spheres with a size distribution in M. d . (Obviously, for sphere packings, particle rotations are 
irrelevant.) Then we describe the solution of the ASC optimization problem using sequential 
linear programming (SLP) techniques. 



A. Adaptive Shrinking Cell (ASC) Scheme for Spheres 



Consider a periodic packing of N spheres with diameters Di, D 2 , . . . , in a fundamental 
cell F of a lattice A in W 1 . As noted in Sec. [Til the fundamental cell is specified by the 
generator matrix Ma = {Ai, A2, . . . , A^}, where Aj are the lattice basis vectors. The design 
variables of the ASC scheme includes both particle displacements, and deformation and 
volume change of the fundamental cell F. Let the global position of the center of sphere 
n (referred to some arbitrary set of non-orthogonal basis vectors) be denoted by x^. For 
simplicity, it is useful to choose orthogonal bases vectors or a Cartesian system. Denote 
by Ax^ a translational displacement of this particle center. Thus, the new position of the 
center of particle n is given by 

K = < + (6) 
However, since the fundamental cell is adaptive, we also need to consider the positions of 
the particle centers in terms of the lattice basis vectors Aj (z = 1, 2, . . . , d). For the center of 
mrticle n, this relative position x^ is related to the global position x^ by the transformation 



44] 



Ma- 1 ■ x* 



(7) 



The adaptive fundamental cell allows for a small change in the fundamental cell 

AM A = £ ■ M A , (8) 
including both volume and shape changes, where e is a symmetric strain tensor, i.e., 



en e i2 ... t\d 

621 622 • • • ^2d 
Ql £d2 ■ ■ ■ £dd 



(9) 
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The new lattice is specified by the new matrix generator 



M A = M A + AM A . (10) 
Substituting the above equation into Eq. (J7|) yields 

x» = M A .x£ = x*+AM A .xi. (11) 

Straining the fundamental cell corresponds to non-trivial collective motions of the particle 
centers. In general, the translational motions of the particles contain contributions from a 
direct part [given by Eq. d^J)] and the collective motion imposed by the adaptive fundamental 
cell. It is this collective motion that enables the algorithm to explore the configuration space 
more efficiently and to produce highly dense packings. 

The displacement pointing from the centroid of sphere m to that of sphere n in terms 
lattice basis vectors \ is given by 

r mn = X n ~~ X m> (12) 

and the counterpart global displacement vector r 9 mn is 

V fnn = X n — X m" (13) 

The Euclidean distance r g mn between the centroids is then given by 

r mn = \ T mn\ = \/ r mn ■ I"mn = \/ r mn ' G ■ Y mn , (14) 

where G = ■ M A is the Gram matrix of the lattice A. Thus, the general mathematical 
formulation of the our ASC optimization scheme for hard sphere packings in d- dimensional 
Euclidean space R d is: 

minimize - 0(x^, . . . , x^; M A ) 

— / \ ( > 

subject to: r 9 mn > D mn , for all neighbor pairs (m,n) of interest. 

where D mn = (D m + D n )/2. 

It should be emphasized that the neighbor pairs here do not necessarily mean nearest 
neighbors. Instead, they are determined by a distance 7 mn , i-e., two spheres m and n are 
neighbors of one another if their pair distance r 9 mn < 7 mn . Thus, by near neighbors of a 
given sphere we mean all of the spheres within some radius of that given sphere. Here, we 
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choose 7 mn = aD mn and a is a positive in the interval [1, L/2D mn ], where L is the length of 
the shortest lattice vector associated with the fundamental cell. For monodisperse packings, 
7 mn is identical for all particle pairs (m, n). For polydisperse packings, 7 mn is generally 
different for each pair (m, n). We call 7 mn the radius of influence sphere associated with 
each pair (m, n). As we will see, these radii are crucial in determining the density and degree 
of disorder of the final jammed states. 



B. Solving the ASC formulation using Sequential Linear Programming 

The formulated problem for the ASC scheme can be solved by considering an equivalent 
sequential linear programming (SLP) problem. Suppose that we only allow the fundamental 
cell to change by a small amount from its original size and shape, and only allow small 
particle displacements; then both the impenetrability constraints and the objective function 
can be linearized. Due to the linearization, one needs now to explicitly places bounds on 
the strain components and the particle displacements, which do not appear in the original 
problem. Once the linearized problem is solved via the linear programming method, a new 
many-particle configuration can be generated using the resulting particle displacements and 
fundamental cell. Then this new configuration is used as an initial configuration, based 
on which a new LP problem is formulated and solved. This sequential linear programming 
procedure is repeated until the increase of the packing density is smaller than some prescribed 
small tolerance value, implying the system is jammed up to high numerical accuracy. 

Consider an initial packing configuration of N spheres within a fundamental cell of volume 
Vp. According to relation (4), the density of the packing is given by 



(16) 



where Vi(Ri) is the volume of sphere % with radius Ri. Now consider a small change of the 
fundamental cell that leads to a cell volume change Av, as well as a change of the positions 
of the centers of every sphere that obey the impenetrability constraint. The density of the 
new configuration in the new fundamental cell is given by 

E,^iW (17) 

If Av is sufficiently small, we can expand as a Taylor series in Aw, keeping only the linear 
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terms, i.e., 



v° F + Av 
v° F 



Av 

F 



= o [l-*r(e)], (18) 

where we have used the relation Av/v° F = tr(e), i.e., the relative small volume change is 
given by the trace of the strain tensor. We see from the above equation that maximizing 
the packing density (ft is equivalent to minimizing the trace of the strain tensor tr(e). 

The nonoverlapping conditions can be linearized in a similar way, which leads to the 
following sequential linear programming formulation for the ASC procedure: 

minimize tr{e) = en + e 2 2 H h 

subject to: 

A . . A/T. . v.A i A ^A _ ci . *-A i . CI . ^ 



Ma • r^ m • e • M A • r* m + Ax* • G • r£ m + Ax* • G • r 



mil 



>l(D 2 mn -^ nm -G-^J+TZ, (19) 
for all neighbor pairs (m, n) of interest, 
Ax% lawer < Ax^ < Ax.% upper , for all n = (1, N), 

flower ^ g <^ g.upper 

Note that the tensor/vector inequalities in (Tl9|) apply to the corresponding components and 
r nm = ~ r mn = X m — X n * s the displacement vector between spheres m and n in the initial 
configuration. The scalar quantity 1Z > is a relaxation variable, which accounts for the 
effects of higher-order terms that would be ignored in the purely linearized problem. The 
value of 1Z is determined by the bounds on the design variables, i.e., 

Tl = max M {-AM. K - Ax^J. (20) 

Introducing TZ has several advantages. Firstly, it makes the linearized problem rigorously 
equivalent to the original (unlinearized) optimization problem. Secondly, it enables us to 
set practically large bounds on the design variables (e.g., the strain components and the 
particle displacements), which in turn enables the algorithm to explore a larger region of the 
configuration space around the initial point and makes it more efficient to generate dense 
packings. In practice, however, we find that as long as the bounds on the design variables 



12 



are chosen to be sufficiently small, one can safely set 1Z = and check the generated packing 
to ensure that no impenetrability constraints are violated. If some of them are violated due 
to the deformable fundamental cell, i.e., leading to a non-positive-definite quadratic term 
associated with the strain tensor in these conditions, the bound widths are reduced to one 
half of their original value and the linearized problem is re-solved. 

Since the objective function and constraints in the new problem are linear functions 
of the design variables, the new problem (fT9l) is solved using standard linear programming 
methods (e.g., simplex method, interior-point method, etc.). Starting from an initial packing 
configuration, we solve the linearized problem and find a new packing, which is denser and 
slightly different from the starting configuration. Then the new packing is used as the 
starting configuration, and a new linearized problem is solved. By repeating this process, 
one actually solves a sequence of linear programming problems to generate a denser packing 
from the previous configuration. 

The choice of the bounds for the design variables is also important in practice. For exam- 
ple, different bound widths for the principle strain components and shear strain components 
correspond to different compression and deformation rates for the packing. Choosing the 
bound widths carefully can dramatically improve the ability of the algorithm to generate 
dense packings. We also emphasize that the choice of neighbor pairs (i.e., the value of 
the influence sphere radius 7 mn ) is also nontrivial, which we will elucidate in the ensuing 
discussion. 



C. Ensuring Jamming 

Finally, we note that our SLP solution procedure guarantees (up to numerical precision) 
that the final locally or globally maximally dense states are indeed strictly jammed under 
periodic boundary conditions. In general, each impenetrability (nonoverlapping) condition 
defines a curved hypersurface in the configuration space that separates accessible and in- 



accessible regions [14|, |45j. Linearization of the impenetrability conditions corresponds to 
replacing the curved hypersurfaces with hyperplanes, which further reduces the accessible 
region (i.e., the linearized conditions are even stronger than the original conditions). How- 
ever, near the jamming point, the available configuration space for the spheres asymptoti- 
cally approaches a closed convex polytope which is exactly determined from the linearized 
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impenetrability conditions |45| . In other words, toward the jamming point, the lin- 
earized conditions become exact nonoverlapping conditions asymptotically. Therefore, the 
SLP algorithm guarantees jamming of the final packing configurations. 

In addition, Our SLP algorithm near the final steps of the process is intimately related 
to a linear-programming (LP) protocol to test for jamming for hard-sphere packings that 
we previously developed 46j. In the LP jamming-test protocol, the interparticle gaps are 
maximized subject to the nonoverlapping conditions. If large gaps can be opened, the 
packing under consideration is not jammed; otherwise it is jammed. For our SLP algorithm, 
the packing density is driven to a local or global maximum, which is not only only collectively 
jammed (since as mentioned earlier, the descent to maxima necessarily involve collective 
configurational motions) but strictly jammed, since the fundamental box is deforming and 
shrinking (on average). In general, if a packing is not jammed, there exist collective particle 
motions and boundary deformations that can lead to a higher density. All of the sphere 
packings in M 2 and M 3 produced by our SLP algorithm are tested using the LP jamming- 
test protocol 46]. The packings with a deformable/shrinking fundamental cell were found to 
be strictly jammed and those with an isotropically shrinking fundamental cell were found to 
be collectively jammed, which verifies the robustness of our algorithm in producing jammed 
packings. 



IV. ENERGY LANDSCAPE, INHERENT STRUCTURES (LOCAL MINIMA) 
AND GLOBAL MINIMA 

A. Energy Landscape Picture 

From a statistical-mechanical point of view, the solution of the ASC formulation using 
a SLP amounts to searching nearby local minima of the energy landscape. The energy 
landscape, defined by the energy (objective function), equal to the negative of the trace of 
the strain tensor associated with the fundamental cell F of the packing or, equivalently, 
the negative of the density 0, is a surface embedded in a (dN + l)-dimensional space, 
where d is the number of degrees of freedom of a single particle and N is the number of 
particles in F. A given initial configuration has an energy that corresponds to a point on 
the landscape. Starting from this point in this multi-dimensional space, our SLP algorithm 
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finds a nearby point, which can be reached from the initial point by a single linear multi- 
dimensional displacement with a lower energy (higher density) subjected to the limits on 
the size and directions of the displacement as well as the impenetrability conditions. This 
procedure is repeated until a jammed state with a high density is reached, which can either 
be a local energy minimum (local density maximum) or a global energy minimum (global 
density maximum) depending on the bound widths on the design variables, as well as the 
influence sphere radii 7 mn for near-neighbor determination. 




(a) (b) 

FIG. 2: (Color online). Influence spheres (i.e., effective interaction range for the central particle in 
so far as the nonoverlapping condition is concerned) and the directions of motion of the particle, 
(a) Only the nearest neighbor of the central (red) particle is included in the influence sphere. The 
central sphere will move away from its neighbor along a line connecting their centers, which allows 
a maximal shrinkage of the fundamental cell. Such a shrinkage would lead to a local maximum 
in the packing density, (b) More near neighbors are included in the influence sphere. The central 
particle moves along a direction that maximizes its distance from all of the neighbors within the 
influence sphere. In other words, the movement of the central particle is coupled with all of the 
near neighbors within the influence sphere and allows a maximal shrinkage of the fundamental cell, 
which after a sufficient number of iterations leads to a global maximum of the packing density. 

As indicated in Sec. III. A, for monodisperse sphere packings, the radius 7 mn of the influ- 
ence sphere for pair (m.n) is identical for each pair of particles. We call this the influence 
sphere for each particle (see Fig. [2]) because only the particles whose centers are within 
the influence sphere will affect the central particle as far as the nonoverlapping conditions 
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are concerned. For polydisperse packings, 7„ m is generally different for each particle pair. 
However, loosely speaking, one can still imagine an effective inffuence sphere around each 
particle, which determines the neighbors that affect the central particle. If the inffuence 
sphere (i.e., 7 mn ) is sufficiently small such that only pairs with minimal separation distance 
are considered, particle motions are controlled by the local environment (see Fig.[2£a)), i.e., 
they move in opposite directions along the line connecting their centroids. If only particle 
pairs (or triplets, etc.) with the minimal separation distance are considered at each stage 
of the SLP procedure, the system rigorously follows the steepest descent trajectory and 
evolves to the associated inherent structure. As the inffuence sphere ( r y mn ) becomes larger, 
successively larger numbers of spheres must respect the nonoverlap condition, all of which 
affect the motion of the central particle (see Fig. [2]^b)). 

Consequently, the inffuence sphere radius 7 mn can be considered to be an effective in- 
teraction range. For each near-neighbor pair, a nonoverlapping condition for that pair is 
included in the SLP formulation, which in general will lead to a different jammed-packing 
solution from one in which that pair of particles are not considered neighbors of one another. 
Therefore, the neighbor pairs effectively interact with each other through the nonoverlap- 
ping condition. A larger 7 mn corresponds to a longer interaction range, which enables more 
directions in the energy landscape to be explored beyond the one that leads to a local den- 
sity increase (i.e., a steepest descent in the energy landscape). We will see in the ensuing 
discussion that this feature of the SLP procedure enables it to find both local and global 
energy minima (density maxima). We note that the number of neighbor pairs n pa i r is a 
monotonically increasing function of 7 mn , namely, n pa i r ~ Imn f° r an y jammed packings 
in M. d . Although a large value of 7 mn is necessary to produce maximally- dense packings, 
one generally only needs an associated small number of particles in the fundamental cell. 
Therefore, the total computational cost when a large value of 7 mn is employed to generate 
maximally- dense packings is usually much smaller than that associated with a small value 
of 7mn to produce disordered packings with a large number of particles per fundamental 
cell. Empirically, we find that when ^ mn > 4D mn , the SLP algorithm generally leads to 
the maximally-dense packings, where D mn = (D m + D n )/2 is the average diameter of two 
spheres. 

Moreover, one can view this SLP compression process physically as a compression of a 
hard-sphere system in a superviscous liquid subjected to periodic boundary conditions in 



1(3 



the absence of a gravitational field. Although the immersion of the hard spheres in the 
superviscous liquid precludes any collisions between the particles, there is a dynamical in- 
teraction with the system boundary. Thus, the only "dynamical" parameters of this physical 
compression are the deformation parameters of the fundamental cell. For this SLP "compres- 
sion" process, the system moves efficiently by collective particle motions to nearby denser 
configurations. 



B. Inherent Structures - Mechanically Stable Local Minima 

When the influence sphere radius 7 mn is sufficiently small (e.g., •j mn ~ 1.5D mn ) and the 
bound widths are taken to be sufficiently large (e.g., |e^-| ~ 0.1 and |Ax A | ~ 0.5D, where 
D is the averaged diameter of all the spheres in the packing), the SLP solution of the ASC 
scheme from some initial configurations leads to a mechanically stable local energy minimum 
(local density maximum), which in principle is the inherent structure associated with the 



starting initial many-particle configuration 
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; see Fig. [3J All initial configurations 



that compress to the same jammed structure, excluding distinctions between states that 
differ by the interchange of identical particles, belong to the same inherent structure. In two 
and three dimensions, the number of distinct inherent structures scales with the number 
of particles iV in the fundamental cell, i.e., exp(aiV), where a is dimensional-dependent 
constant of order unity ^7|. 

Stillinger, DiMarzio, and Kornegay were the first to introduce the idea of inherent struc- 
tures for two-dimensional hard circular disks 49j. However, at that time they did not 
explicitly use the terminology "inherent structure." Nonetheless, they indeed proposed a 
conceptual procedure to produce inherent structures for hard-sphere systems and applied 
it to hard disks in two dimensions. They considered a dilute configuration of N spheres 
centered at x 1; . . . , Xjv, and subsequently expanded the spheres. During such an expansion, 
there will occur a first contact between the pair of spheres with the minimal original pair 
separation distance. As the spheres are further expanded, one should move this pair apart 
(each member at an equal rate along the line joining the centers) just to maintain the con- 
tact. Other pairs in the system similarly will come into contact, and will be rearranged by 
the same prescription. As the process proceeds, triplets, and larger sets of spheres will touch, 
and subsequently should be moved by a generalization of the pair procedure. These authors 
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FIG. 3: (Color online). A schematic diagram illustrating the idea of inherent structures for a system 
composed of N hard spheres as adapted from Ref. 49(, though the term "inherent structure" was 
not explicitly used there. The horizontal axis labeled ~K N stands for the entire set of centroid 
positions, and the density <j) increases downward. The jagged curve is the boundary between 
accessible (upper) configurations (shown in yellow or light gray in the print version) and inaccessible 
(lower) configurations (shown in blue or dark gray in the print version). The deepest point of the 
accessible configurations corresponds to the maximal density packings of hard spheres. 

proposed that this procedure can be mapped into a steepest-descent problem. Namely, if one 
has a set of n spheres in contact at xi, . . . , x n , the displacements dxi, . . . , <ix n , which main- 
tain contacts under the diameter increase AD, should be selected to minimize the positive 
definite form 

JL /A Y .\ 2 

(21) 



i=l 



subject to those constraints. In particular, these authors commented that "for very large 
N, the chance of selecting (at random) a set of initial positions that would lead to jamming 
in the regular close-packed array becomes very small." In modern language, this statement 
means that random initial configurations for the case of identical spheres would lead to 
maximally random jammed states for d > 3, which are isostatic configurations {50I, with 
very high probability. However, for two-dimensional monodisperse circular disk packings, it 
is well known that random initial configurations lead to highly ordered crystalline packings 
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with very high probability jl^] . If an entropic measure of disorder was employed, such two- 
dimensional crystalline packings would be designated to be the most disordered packings or 
MR J states, which is clearly incorrect. This is one of the many reasons that an entropic (i.e ., 
occurrence frequency of the same configurations) measure of disorder can be misleading 14]. 

It is not difficult to see that our SLP algorithm when provided with the appropriate 
parameters closely follows the aforementioned conceptual steepest-descent procedure for ob- 
taining inherent structures. Specifically, if the non-overlapping conditions are only specified 
for nearest-neighbor pairs, triplets etc., i.e., a small influence sphere radius 7 mn is used (see 
Fig. [2]), the SLP algorithm will only move apart those spheres that are explicitly considered 
in the non-overlapping conditions such that the largest possible compression of the system 
(consistent with the non-overlapping conditions and the bounds) can be performed (through 
applying the solved strain tensor). In other words, starting from an initial configuration, the 
SLP algorithm can lead to a maximal density increase subject to the constraints imposed 
by closest neighbors, following the trajectory of steepest descent in the energy landscape. 
This is an equivalent re-formulation of the Stillinger-DiMarzio-Kornegay steepest-descent 
procedure. Importantly, the SLP algorithm is not limited to steepest-descent mappings. We 
will see in the next section that by using relatively large influence spheres, one can access 
unusual inherent structures from random initial configurations. 



We note that Zinchenko 



5l| proposed and implemented an numerical algorithm concep- 



tually equivalent to the Stillinger-DiMarzio-Kornegay steepest descent procedure to produce 
disordered jammed spheres packings in three dimensions. However, the Zinchenko algorithm 
involves the numerical solution of many different differential equations (dN of them) as well 
as the solutions to dN algebraic equations that prescribe the densification process that al- 
lows the particles to swell while retaining sphere contacts to the extent possible. This makes 
the algorithm both computationally expensive and incapable of generating ordered pack- 
ings. Furthermore, O'Hern et al. employed a conjugate-gradient (CG) method to produce 
jammed packings of soft spheres interacting with short-ranged repulsive power-law pair po- 



tentials 52J. Starting from a dilute initial configuration, the spheres are allowed to grow, 



which cause an increase of the total energy of the system. Then CG method is used to 
relax the system to zero energy. The procedure is repeated until the total energy cannot be 
relaxed to zero, and the system is considered jammed. Although not explicitly indicated, 
this algorithm in general produces inherent structures of such soft-sphere systems. 
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C. Energy Minima - Density Maxima 



Though our SLP algorithm can produce inherent structures, it is crucial to emphasize that 
our algorithm is much more general than the aforementioned conceptual steepest-descent 



procedure [49|, the Zinchenko protocol [51] and the conjugate- gradient method [52|. In 
particular, including more neighbors associated with any given sphere in the set of non- 
overlapping conditions, i.e., increasing the radius of the influence sphere 7 mn , effectively 
introduces "long-range interactions" (i.e., 7 mn ~ 3.5D mn ). Together with smaller bound 
widths (e.g., | e^- 1 ~ 0.01 and |Ax A | ~ 0.05-D (where D is the effective diameter of all 
of the spheres in the packing), the SLP solution of the ASC scheme from random initial 
configuration can lead to unusual inherent structures, such as a mechanically stable global 
energy minimum (global density maximum). 

Recall that when the radius of influence sphere (j mn ) becomes larger, successively larger 
numbers of particles must respect the nonoverlap condition, all of which affect the motion of 
the central particle. A sufficiently large 7 mn implies that the particles experience effectively 
long-range "interactions" via the influence spheres, which enable more directions in the 
energy landscape to be explored besides the steepest-descent direction. In other words, 
multiple minima can now be explored; and instead of moving to a local energy minimum 
(local density maximum) by a steepest descent trajectory, the system is able to evolve toward 
the deepest minimum available until it is in the basin of the global minimum. 

D. Unusual Inherent Structures - Low-Density Jammed Packings 

From its definition, it is clear that a inherent structure is highly dependent on its asso- 
ciated initial configurations. Therefore, an innovative procedure to generate unusual initial 
configurations in principle enables one to obtain unusual inherent structures. Moreover, 
there is no reason to limit oneself to random initial configurations alone. For example, in 
Ref. 4l| diluted MRJ packings are employed as initial configurations to produce low-density 
jammed packings in R 3 . 
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V. JAMMED DISORDERED AND ORDERED PACKINGS 



In this section, we employ our SLP algorithm to produce maximally random jammed 
(MRJ) packings as well as maximally dense packings of hard spheres in <i-dimensional Eu- 
clidean space M. d for d = 2, 3, 4, 5 and 6. In M 1 , there is only a single jammed state, namely, 
the integer-lattice packing with unit density. Therefore, we will not consider this trivial case 
here. The characteristics of the packings produced by the SLP algorithm are compared to 
those obtained using the LS algorithm, which verifies the robustness and low-computational 
cost of our SLP algorithm. Moreover, we show that the SLP algorithm is superior to the LS 
algorithm in several aspects, especially in producing maximal dense packings. 

A. Maximally Random Jammed Packings 



FIG. 4: (Color online). Maximally random jammed packings in low dimensions via SLP algorithm, 
(a) A MRJ packing of 1000 monodisperse spheres in three dimensions, (b) A MRJ packing of 500 
polydisperse circular disks in two dimensions. 

Random dilute packings with a density <p — 0.05 are used as initial configurations to 
produce MRJ packings for d — 2, 3, 4, 5 and 6 using the SLP algorithm. In two dimensions, 
polydisperse disks with diameters uniformly distributed within [D, 2D] are used to generate 
MRJ packings because monodisperse circular disks have an inevitable propensity to form 




(a) 



(b) 
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FIG. 5: (Color online). Pair correlation function g<i of maximally random jammed packings of 
monodisperse hyperspheres in selected dimensions via the SLP and LS algorithm. 



highly ordered packings in the jamming limit, as discussed in Sec. IIVBI In higher dimensions, 
monodisperse spheres are used. As indicated in the previous section, to produce MRJ 
packings (which are the maximally disordered inherent structures) , a small influence sphere 
radius 7 mn is used such that evolution of the packing in the energy landscape follows a path 
of steepest descent (e.g., ■y mn ~ 1.5D mn , \eij\ ~ 0.1, and |Ax A | ~ 0.5D). The number of 
particles in the fundamental cell are N = 500,2000,4000,8000, 12000 respectively for d = 
2, 3, 4, 5, 6. The packing is considered jammed and the simulation is terminated if the increase 
of the packing density is less than 10~ 8 . Figure H] shows typical packing configurations in 
two and three dimensions. 

The characteristics of the monodisperse packings are compared to those obtained using 
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TABLE I: Characteristics of MRJ packings produced using the SLP and LS algorithms. <j) is the 
density of final jammed packing and f r is the fraction of the rattlers. t s is the total simulation 
time. For LS algorithm, t s includes the time of simulations with both the initial large expansion 
rate and the fine-tuning expansion rate. 



LS algorithm 



SLP algorithm 



d=3 = 0.642 ± 0.005 f r = 0.030 ± 0.003 

t s = 1.5 hrs 
d=4 (j) = 0.460 ± 0.006 f r = 0.012 ± 0.002 

t 8 = 4.8 hrs 
d=5 (j) = 0.315 ± 0.005 f r = 0.008 ± 0.002 

t 8 = 14 hrs 

d=6 (j) = 0.201 ± 0.005 f r = 0.005 ± 0.001 
t s = 193.5 hrs 



0.640 ± 0.004 f r = 0.028 ± 0.003 

t s = 10 mins 
0.465 ± 0.005 f r = 0.014 ± 0.002 

t s = 46 mins 
0.310 ± 0.005 f r = 0.006 ± 0.002 

t s = 3.2 hrs 
0.199 ± 0.005 f r = 0.005 ± 0.001 

t s = 8.3 hrs 



the LS algorithm 25], which has been verified to be a robust protocol to produce MRJ 
packings. A two-step procedure is used for the LS algorithm to generate jammed disordered 
packings: a large initial expansion rate of ~ 0.01 and a small fine-tuning expansion rate 
of 10~( 3+d ) (where d is the spatial dimension). A comparison of the packing density and 
fraction of rattlers (i.e., movable particles caged by their jamming neighbors) is given in 
Table I. Each density <fi and the fraction of rattlers f r are obtained by averaging over 5 
configurations. We note that the MRJ densities obtained via the LS and SLP algorithms 
are very close to each other. We expect that in the infinite-system limit, the MRJ densities 
generated via the two algorithms should be essentially identical to one another. Figure [5] 
shows the pair correlation functions gi of the packings obtained using both the SLP and LS 



algorithm 



55[ . It can be seen clearly that the MRJ packings produced by our SLP algorithm 



are consistent with those produced by the LS protocol across dimensions. Specifically, in high 
dimensions, the split-second peak in <?2, present for d = 3, gets dramatically diminished and 
oscillations in g% get significantly dampened. These findings are consistent with a recently 
proposed "decorrelation principle" 2^] that states that unconstrained spatial correlations in 
disordered packings should vanish asymptotically in the high-dimensional limit. The two- 
dimensional polydisperse packing possesses a density <p = 0.846 and approximately 3% of 
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rattlers. 

We note that in producing MRJ packings using the LS algorithm, a large expansion rate 
is necessary in the beginning of the simulation to drive the system out of equilibrium but 
it is undesirable toward the end, which may lead to unjammed configurations. Therefore, 
a variable expansion rate is needed and must be very small at the end of the simulation 
to ensure jamming, as we have done here. In our SLP algorithm, the control parameters 
are single- valued. Moreover, the SLP algorithm naturally leads to jammed packings, as we 
discussed in Sec. III. B, which makes it superior to the LS algorithm in producing jammed 
packings. 



B. Maximal Density Packings 

In order to further test the SLP algorithm, we apply to see whether it can recover the 
densest packings of congruent spheres in two and three dimensions and densest known 
congruent sphere packings in R 4 , R 5 and R 6 . Recall that we do not consider the case in R 1 , 
because that corresponds to the trivial solution of the integer-lattice packing. In R 2 and 



d3 



the densest pac 



dngs of congruent spheres have been proved to be the triangular lattice 



[531 ] and fee lattice [32| , respectively. One of the generalizations of the three-dimensional fee 
lattice to higher dimensions is the d- dimensional checkerboard lattice D^, defined by taking a 
hypercubic lattice and pla cing spheres on every site at which the sum of the lattice indices is 
even, i.e., every other site 22J. The four- and five-dimensional checkerboard lattices -D4 and 
D§ with densities (p = 7r/16 = 0.6168 . . . and <\> = ^/2ir 2 /30 = 0.4652 . . ., respectively, are 
believed to be the densest sphere packings in those dimensions. For d = 6, the densest sphere 
oacking is conjectured to be the "root" lattice Eq with density = "v/3tt 3 /144 _ q.3729 . . . 



22]. 



It is important to note that in order to generate perfect crystalline jammed packings using 
a finite simulation cell with periodic boundary conditions, the number of particles in the 
fundamental cell must be consistent with its shape. For example, if a fixed cubic fundamental 
cell or rhombohedral fundamental cell (associated with the fee lattice) in three dimensions 
is used, this "magic" number of spheres should be N = 4m 3 and N = m 3 (m = 1, 2, 3 . . .), 
respectively. In most cases, the fundamental cell associated with the densest packing is 
not known a priori and hence one needs to explore a variety of different particle numbers 
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TABLE II: Characteristics of densest jammed packings produced using the SLP and LS algorithms. 
Here is the density of final jammed packing, N is the number of particles in the fundamental 
cell, and £, is the total simulation time. 



LS algorithm 



SLP algorithm 



d=2 



N 



= 0.9068 . . 
0.1 mins 
d=3 N = 27 = 0.7408 . 

t s = 6.5 mins 
d=4 N = 81 = 0.5608. 

t s = 1.6 hrs 
d=5 iV = 243 = 0.4154. 

t s = 47.5 hrs 
d=6 iV = 729 = 0.3287. 

t a = 283.5 hrs 



N = 9 

t a -~ 

N = 27 
N = 81 

t a = 

N = 243 

ts = 

N = 729 



j> = 0.9068. . . 
- 0.1 min 
(j) = 0.7408. . . 
: 1.5 min 
= 0.6168... 
4.5 mins 
= 0.4652.. 
: 11 mins 
= 0.3729.. 



t s = 27 mins 



(by carrying out many simulations) to obtain the maximal density packing with typical 
protocols that use a fixed cell. Therefore, it is highly desirable to explore systematically the 
shape of the appropriate fundamental cell by employing a single simulation with variable 
particle numbers, which our SLP algorithm enables one to do. In particular, because our 
SLP algorithm does not simulate the real many-particle dynamics, a very small value of 
iV can be used, which can significantly reduce the simulation time. For example, one can 
choose iV to be as small as one, which reduces the SLP algorithm to an lattice-optimization 
algorithm. For the LS protocol, using a very small value of N will introduce systematic 
errors due to lack of collisions between the particles within the fundamental cell or because 
collisions are not properly treated. 

To produce maximal-density packings in M d , initial packing configurations with N = 3 d 
spheres in (hyper) cubic fundamental cells possessing a density lower than the "equilibrium" 
hard-sphere freezing density (Ref. |26j) for each dimension (i.e., = 0.45,0.4,0.3,0.15,0.06 
for d = 2,3,4,5,6, respectively) are used for both the LS and SLP algorithm. We note 
that for the current implementation of the LS algorithm, a deformable fundamental cell 



program is available for d = 2 and 3 36]; for d > 4, only a fixed cubic fundamental cell 
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is available 



26(. The sphere expansion rate used to generate a maximally dense packing 



in R d is 10 ( 3+d ). For the SLP algorithm, a deformable fundamental cell is used in all 
dimensions. The influence sphere radius is chosen to be 7 mn ~ 3.5D), and |ejj| ~ 0.01 
and |Ax A | ~ 0.05D are used. The final densities produced by the two algorithms and the 
corresponding computational times to produce the jammed states are compared in Table II. 
It can be seen that the SLP algorithm successfully reproduces the densest known packings 
in all dimensions; while our implementation of the LS algorithm only was able to reproduce 
the maximal-density packings in R 2 and IR 3 . Although the LS algorithm does not generate 
perfect crystals for d = 4, 5 and 6, the high densities, it produces suggest that these packings 
contain very large crystallized regions. We have also used iV = 128 in M 4 and N = 512 in 
IR 5 for the LS algorithm, which are the magic numbers for cubic fundamental cells in four 
and five dimensions, respectively, as well as an expansion rate of 10 -10 . In M 4 , the packing 
converges to the ^-lattice packing with (ft = 0.6168. . .; while in M 5 , the final packing is 
partially crystallized with (ft = 0.4392 .... 

These results clearly demonstrate the advantage of our SLP algorithm in producing max- 
imally dense sphere packings, i.e., the SLP algorithm is both much more computationally 
efficient and more robust in finding the densest known sphere packings than the LS algo- 
rithm. In particular, the adaptive fundamental cell is crucial for obtaining perfect maximal 
density packing structures. In addition, we have tried a variety of particle numbers in the 
fundamental cell (e.g., 1 < N < 2 d for R d ) in four, five and six dimensions and have not 
obtained any packings possessing a larger density than the densest known packings. This 
provides numerical support that the densest known packing in these dimensions may indeed 
be the conjectured maximal- density packings. 



We note that very recently Kallus, Elser and Gravel 



54J applied a periodic "divide" and 



y dense packings of hard particles. In this 



"concur" (PDC) algorithm to generate maximal . 
variant of the PDC algorithm, they follow Ref. 43J by allowing the fundamental cell to de- 
form and change shape during the simulation. Interparticle impenetrability constraints are 
relaxed in the beginning of the simulation such that overlapping intermediate configurations 
are allowed. Towards the end of the simulation, all overlaps are systematically removed. 
Among other applications, these authors applied the PDC algorithm to successfully repro- 
duce the densest known lattice sphere packings (not necessarily the densest packings) and 
the best known lattice kissing arrangements in up to 14 and 11 dimensions, respectively, 
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since they were limited to using only a small number of particles per fundamental cell. Al- 
though the PDC algorithm is efficient in producing dense lattice sphere packings, it does 
not appear to be capable of generating jammed sphere packings with a diversity of disorder 
and density that can be produced via either the SLP or LS algorithm. 

VI. THE ORDER MAP 




0i 1 , 1 , 1 1 

0.5 0.6 0.7 

Packing Density 

FIG. 6: Order map for strictly jammed packings of congruent spheres in three dimensions with 
N ~ 1000 generated using the SLP algorithm: translational order metric T versus packing density 
cf). The black circles are data produced via the SLP algorithm using random initial conditions. 
While the density in this data set ranges essentially continuously from <f> ~ 0.64, which coincides 
with that of the MR J state, up to the maximal density of 4> = 7r/\/l8 = 0.7408 . . ., the degree of 
order at fixed density is variable. The red squares are data corresponding to disordered jammed 
packings with anomalously low densities, with density as low as <f> ~ 0.595, which were generated via 
the SLP protocol using special initial conditions involving diluted MRJ packings, as described in the 
text. The dashed blue curves represent the spectrum of packings generated by randomly filling the 
vacancies in three of the tunneled crystal packings down to the density (f> = V2ir/9 = 0.49365 . . . 
(green triangles), each of which leads to the corresponding maximal-density Barlow packing [8] 
(with cp = 71-/^/18 = 0.7408 . . . ) 
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Besides the capability of generating MR J packings and maximally dense jammed packings 
efficiently and robustly, the SLP algorithm has the capacity to produce jammed sphere 
packings of with a wide range of densities and degrees of disorder/order, i.e., a diverse set 
of inherent structures. Here we illustrate this capability by focusing on identical spheres 
in three dimensions. By tuning the bound widths and the influence sphere radius 7 mn , 
the density at which the packing is jammed and its degree of disorder can be controlled 
to a great degree. In particular, the ranges of the parameter values used are as follows: 
7 mn G (1.5D mn ,3.5D mn ), \e i:j \ G (0.01,0.1) and |Ax A | G (0.05D, 0.5D) (where D is the 
diameter of the spheres in the packing). In general, a larger value of 7 mn and smaller bound 
widths allow for a larger number of particles to be effectively coupled to one another, which 
in turn leads to a higher jammed packing density. By continuously varying the values of 
the control parameters, a spectrum of packing density <p £ [0.64, 0.74048...] can be obtained 
when random dilute packings are used as initial configurations. This density range is also 
achievable by the LS algorithm, but we shall see that disordered jammed packings with 
densities significantly lower than 0.64 can be generated using the SLP algorithm. 

Typical configurations of such packings with G [0.64, 0.74048...] and iV »s 1000 are 
mapped onto a density-order-metric diagram (black circles). The diagram or order map em- 
phasizes a "geometric-structure" approach to analyze pac king s by characterizing individual 
configurations, regardless of their occurrence probability 



results 



14J. Here the translational order 



metric % 4|is used to quantify the order of the pac 
3l|). The translational order metric T 



dngs (other order metrics yield equivalent 
4j is defined as 



T 



N s 



N s 



n FCC 



yjideah 



(22) 



where rii is the average occupation number for the shell % centered at a distance from a 
reference sphere that equals the zth nearest-neighbor separation for the open FCC lattice 
at that density and N s is the total number of shells for the summation; nf eal and nf are 
the corresponding shell occupation numbers for an ideal gas (spatially uncorrelated spheres) 
and the open FCC lattice. For a completely disordered system (e.g., a Poisson distribution 
of points) T = 0, whereas T = 1 for the FCC lattice. 

Moreover, we have produced jammed disordered packings with densities as low as 0.595, 
which is an anomalously low value, using the SLP algorithm from special initial configura- 
tions. Specifically, a small fraction (f s = 0.1% — 2.5%) of particles are removed from MRJ 



28 



packings which are obtained using the SLP algorithm to produce unjammed initial configu- 
rations. The remaining unjammed spheres are then compressed to a jammed state using the 
SLP algorithm. Such a procedure is repeated n r = 5 — 10 times before a lower limit on the 
final jammed packing density <j) is reached. Low-density jammed packing with densities in 
the range G [0.6, 0.64] have been produced (removing the rattlers in the packing results 
in a slightly lower density, <p ~ 0.595). These packings are mapped onto the 0-T diagram 
(red squares). We note that a similar procedure was employed in Ref. 4JJ to produce low- 
density jammed packings. However, there the LS algorithm instead of the SLP algorithm 
was used to generate the initial configurations. Moreover, we would like to emphasize that 
the low-density jammed packings produced by the SLP algorithm are completely unrelated 
to the so-called ill-defined "random loose packings" of hard spheres, which have been shown 
to be not even strictly jammed 46]. 



"Tunneled" crystal packings (green triangles) [8( and maximal density packings (fee, 
hexagonal close-packed, and their stacking variants called the Barlow packings) are also 
shown on the order map. The blue dashed curves show the spectrum of packings generated 
by randomly filling the vacancies in the tunneled crystal packings, which leads to the cor- 



responding maximal density packings. These dashed curves were first reported in Ref. [41] . 
The tunneled crystals have density <p — V2k/9 = 0.49365. . . are currently the best candi- 
dates for the lowest density strictly jammed packings in three dimensions, bf Each sphere 
in any of the tunneled crystals contacts exactly 7 others and therefore such packings are 
hyperstatic. 

We note that for comparison to existing literature, the data presented in Fig. for the 
order map include the rattlers, which constitute a small percent of spheres near the MRJ 
state (e.g., 2 ~ 3%). The fraction of rattlers decreases as the packing density increases or 
decreases from MRJ state along the lower boundary of jammed states (in either direction 
in density). In addition, consistent with the results reported in Ref. 41|, we find that 
the average contact number Z per particle increases from the isostatic value of 6 at the 
MRJ state as the density increases or decreases from the MRJ-state value along the lower 
boundary of jammed states. In particular, at <fi ~ 0.6 the average contact number Z « 6.3. 
It is expected that Z will increase to a contact value per particle of 7 as is decreased to 
the tunneled-crystal value <p = 0.49365 . . .. On the other hand, as <fi is increased from the 
MRJ value, the average contact number Z continuously increases from 6 all the way up to 
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12 for the hyperstatic Barlow packings with a density of 0.74048 .... Thus, moving off the 
MR J state along the lower boundary of jammed states is associated with an increase in the 
average contact number Z and the degree of order. 



VII. DISCUSSION AND CONCLUSIONS 

In this paper, we have proposed and implemented a sequential-linear-programming pro- 
cedure to solve the adaptive-shrinking-cell optimization problem to generate jammed sphere 
packings. The SLP procedure is particularly suitable and natural for sphere packings with 
a size distribution because the objective function and impenetrability constraints can be 
exactly linearized and the final state is ensured to be jammed in principle. We have shown 
that the SLP algorithm can produce robustly a wide spectrum of jammed sphere packings 
in R d for d = 2, 3, 4, 5 and 6 with varying degrees of disordered and densities. In particular, 
we applied the algorithm to generate various disordered packings as well as the maximally 
dense packings for d = 2,4,5 and 6. Moreover, we showed that our SLP algorithm can pro- 
duce with high probability a variety of collectively jammed packings with a packing density 
anywhere in the wide range [0.6,0.7408. . .] in three dimensions, which supports the view 
that there is no universal jamming point that is distinguishable based on the packing den- 
sity alone Q, [31], 41 1. Our jammed sphere packings are characterized and compared to the 
corresponding packings generated by the well-known and versatile Lubachevsky-Stillinger 
molecular-dynamics packing algorithm. Compared to the LS procedure, our SLP protocol is 
able to ensure that the final packings are truly jammed, produces disordered jammed pack- 
ings with anomalously low densities, and is appreciably more robust and computationally 
faster at generating maximally dense packings, especially as the space dimension increases. 
Therefore, the SLP algorithm for producing jammed sphere packings retains the versatility 
and advantages of the LS algorithm while improving upon its imperfections. 

An important feature of the SLP algorithm is that it can produce a broad range of inher- 
ent structures (locally maximally dense and mechanically stable packings), besides the usual 
disordered ones (such as the maximally random jammed state) with very small computa- 
tional cost compared to those of the best known packing algorithms by tuning the radius 
of the influence sphere. This is to be contrasted with other packing protocols, such as the 



conjugate- gradient method for soft three-dimensional spheres 52| and the Zinchenko algo- 
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rithm for three-dimensional hard spheres 51], which can only lead with high probability to 
the MRJ-like inherent structures from random initial configurations. These protocols are 
limited in the sense that the trajectory from one point to another on the energy landscape 
is strongly limited, i.e., only those directions corresponding to local steepest descents can 
be chosen. This makes the final packings strongly dependent on the random initial config- 
urations. As we have shown, exploring other directions in the energy landscape can lead to 
a variety of jammed packing structures with variable disorder and density even if the same 
random initial configuration is used. In addition, it is not clear that these other packing 
protocols would lead to inherent structures with the same level of order/disorder, since no 
order metrics have been explicitly measured for the resulting jammed packings. This em- 



phasizes a major point made in Refs. 



Q,El, 



4l| , namely, although the diversity of packing 



protocols have yet to be fully explored, there is already substantial evidence that available 
algorithms, including the one reported in this paper, can produce jammed packings at the 
same density but with a wide range of order/disorder. 

A straightforward generalization of the current SLP procedure for hard-sphere packings in 
Euclidean space M d is to devise a linear-programming protocol for packing spheres in curved 
spaces, such as the surface of a sphere. In such cases, additional constraints on the sphere 
displacements are necessary to keep them on the spherical surface, but these constraints are 
simple to incorporate. 

One advantage that the LS packing protocol has over any SLP algorithm is in producing 
jammed packings of smoothly shaped nonspherical particles (e.g., ellipsoids and superballs) 
because the nonoverlap functions are nonlinear 36M39I] and linearization of the constraints 
is no longer rigorous and hence jamming cannot be guaranteed. However, in such instances, 
the ASC scheme can be solved using Monte Carlo methods, as was done for hard-polyhedron 
packings 43j, or by using nonlinear programming methods 56[. When the particle shape 
deviates only slightly from a perfect sphere, a LP solution could still be possible. Moreover, 
if one is interested in the true dynamics that leads to jamming, the LS algorithm is clearly 
preferred over the SLP protocol, which is a deterministic algorithm that is incapable of 
capturing the real dynamics. However, if one is simply interested in generating jammed 

packings without regard to history, which the geometric-structure point of view advocates 

i 14, 0], then the SLP protocol is preferable. 



Finally, we note that it is highly desirable to devise packing protocols that targets the 
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density or other packing characteristics while ensuring] amming. It is natural to formulate 



such an inverse problem as an optimization problem [57J and to solve it for sphere packings 
using a LP solution procedure. In future work, we will generalize the current SLP algorithm 
to create even lower density jammed states of hard spheres than have been able to produced 
thus far. In three dimensions, this could enable the large gap without any jammed states 
between densities ~ 0.49 and ~ 0.6 with relatively low order metric values (see Fig. [6]) to 
be filled in. 
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